This chapter shows the capabilities of the simulator and how it was implemented and how it could be used. It presents the algorithms for certain interesting parts of the program so that the reader would be able to recreate it if they wished replicate the experiments performed for this project.

For the system I needed to develop in order to run the experiments, I obtained a list of requirements that it had to achieve:
\begin{itemize}
\item Allow the user to determine the grid size, number of steps and the initial population size, birth probability and death probability of multiple species.
\item Switch between discrete step cellular automata and continuous time Gillespie simulator via flag.
\item Switch between the different discrete step rules as specified in Section \ref{sec:back.how.disc} via flag.
\item Output results in a normalised form (i.e. have the results for the cellular automata be comparable to the Gillespie simulator).
\item Be capable of running from a command line interface.
\item Be efficient in regards to the speed for running the simulations.
\item Have a GUI which can be turned on or off via a flag.
\end{itemize}

\section{Using the simulator}
\label{sec:meth.using}
The system allows for the user to set the rules and the initial state of the simulation via command line arguments. In order to run the simulator, the user must provide the dimension of the grid and the number of steps. For every species (up to a maximum of ten), the user must also specify the initial population, the birth probability and the death probability of that species. For example, to run a simulation with a grid of size 100*100, for 1000 steps with two species, one with an initial population of 60, birth probability of 0.4 and death probability of 0.2, and the other with an initial population of 45 with birth probability of 0.35 and death probability of 0.1, the user would enter "100 1000 60 0.4 0.2 45 0.35 0.1". This would then run the simulation with this setup using the default discrete time step rules (see Section \ref{sec:meth.using.rules}), and output for each step, the step number, the population of each species and then the total population.

There are also optional flags which the user can include, at any place in the arguments, which are:
\begin{itemize}
\item \textbf{`-g'} - Turns on the Grid part of the GUI (see Section \ref{sec:meth.GUI.grid})
\item \textbf{`-c'} - Turns on the Chart part of the GUI (see Section \ref{sec:meth.GUI.chart})
\item \textbf{`-m'} - Counts the number of patches, the largest patch size (see Section \ref{sec:exp.rux.spat}) and Moran's I (see Section \ref{sec:exp.rux.pers}) for every step.
\item \textbf{`-r'} - Allows the user to specify which rules to use when discrete step cellular automaton, or whether to run the Gillespie simulator.
\item \textbf{`-s'} - Allows the user to set the "seed value" for the randomisation function in the program (which will allow for results to be repeated).
\end{itemize}

\subsection{Rules}
\label{sec:meth.using.rules}
The system runs the discrete step cellular automata with the rules set to Sequential, birthFixed and time type 2 by default. The user can run any other combination of rules using the flag `-r' followed by the a string consisting of the shorthand keys for each rule (as shown in Section \ref{sec:back.how.disc}) in the order of Spatial, Temporal and then Time Type. For example, for Random deathConditional with time type 1, the user would use the argument `-r RdC1' and for Sequential bothRandom with time type 2, the user would use `-r SbR2'. To get the system to run as a Gillespie simulator the argument would be `-r gil'.

\section{Implementation}
\label{sec:meth.impl}
The grid is represented using a one-dimensional integer array to hold the position of each organism within the program. Each species is represented by a different number and 0 represents an empty cell. A list of indexes for each of the eight neighbours (the Moore neighbourhood) is calculated for every grid position at the beginning of the program and stored in an array as well. This method was used rather than using a two-dimensional array to represent the grid and modular arithmetic to calculate the neighbourhood as it allows for only one array index to be needed in order to access any cell in the grid and also allows for the neighbourhood of each cell to only have to be calculated once and then only an array index is required to retrieve the information after that. This allowed for the simulator to execute much faster as for an example simulation with a grid size of 100 by 100, a two-dimensional array would require two array indexes for every cell in the grid (20,000 array indexes per step) as well as performing the modular arithmetic for every birth. Using the implementation I chose meant that there would only be one array index per cell in order to scan through the grid, with another index to obtain all of a specific cells neighbours in the event of a birth.

The details for each species (the population size and birth and death probabilities) are also stored in arrays which are indexed using the value used to represent it in the grid. The program uses a pseudo-random number generator so that simulations can be repeated exactly. This way, if a lot of simulations need to be run (or simulations with a large number of steps), the large output files do not need to be stored and instead, if the arguments that the simulation was run with and the seed number are stored, then the simulation can be re-run to obtain the results.

\subsection{Cellular Automata}
\label{sec:meth.impl.ca}

\begin{figure}[H]
\centering \small
\begin{footnotesize}
\begin{tabular}{ll}
A1 & \textbf{for} cell \textbf{in} grid\\
A2 & \hspace{0.25cm} \textbf{if} cell is not empty\\
B1 & \hspace{0.5cm} \textbf{if} randomDouble $\le$ birthProb\\
B2 & \hspace{0.75cm} \textbf{if} randomlySelectedNeighbour is empty:\\
B3 & \hspace{1.0cm} randomlySelectedNeighbour = \textbf{not} empty\\
B4 & \hspace{1.0cm} speciesPopulation = speciesPopulation + 1\\
C1 & \hspace{0.5cm} \textbf{if} randomDouble2 $\le$ deathProb\\
C2 & \hspace{0.75cm} cell = empty\\
C3 & \hspace{0.75cm} speciesPopulation = speciesPopulation - 1\\
\end{tabular}
\end{footnotesize}
\caption{General pseudocode for discrete step cellular automata using the birthConditional rule for the temporal sequence}
\label{fig:disc}
\end{figure}

Figure \ref{fig:disc} shows the pseudocode for each step in the discrete step cellular automata using birthConditional rules. Differences in the ordering for the birth and death evaluations (sections B and C respectively) would be applied depending on the rules set by the user for the simulation. For a Fixed sequential ordering, the loop on line A1 in Figure \ref{fig:disc} would be repeated twice, once for the birth evaluations and once for the death. In order to decide whether an event will occur (lines B1 and C1 in Figure \ref{fig:disc}), a random number between 0 and 1 is chosen and the event only occurs if the random number is less than or equal to the probability of that event. For a birth event, the site for the new organism is chosen by randomly selecting one of the cells within the neighbourhood of the parent and then the birth will only occur if that cell is not currently occupied. To iterate through the grid for each step (line A1 in Figure \ref{fig:disc}), an array is kept containing the indices of the grid in order. If the Random spatial option is set, then the array of indices is iterated through and swapped with a randomly selected number in the array for every step in the simulation i.e.:

\begin{figure}[H]
\centering \small
\begin{footnotesize}
\begin{tabular}{ll}
A1 & \textbf{for} index \textbf{in} array\\
A2 & \hspace{0.25cm} choose random position in array\\
A3 & \hspace{0.25cm} swap the values from the index and the random position in the array\\
\end{tabular}
\end{footnotesize}
\caption{Pseudocode for randomising the spatial order}
\label{fig:rand}
\end{figure}


\subsection{Gillespie Simulator}
\label{sec:meth.impl.gs}

\begin{figure}[H]
\centering \small
\begin{footnotesize}
\begin{tabular}{ll}
A1 & \textbf{for} probabilities of all events\\
A2 & \hspace{0.25cm} averageRate = convertToRate(probability) * populationSize\\
A3 & \hspace{0.25cm} totalRates = totalRates + averageRate\\
B1 & \textbf{for} all averageRates\\
B2 & \hspace{0.25cm} prob = averageRate/totalRates\\
C1 & tau = time until next event\\
C2 & Time = Time + tau\\
C3 & event = randomDouble\\
D1 & \textbf{for} all events\\
D2 & \hspace{0.25cm} \textbf{if} event $\ge$ eventProb$_{i}$ \textbf{and} event $<$ eventProb$_{i+1}$\\
D3 & \hspace{0.5cm} eventProb$_{i}$ occurs\\
\end{tabular}
\end{footnotesize}
\caption{General pseudocode for Gillespie's Algorithm}
\label{fig:gill}
\end{figure}

Figure \ref{fig:gill} shows the pseudocode for Gillespie's algorithm. For this implementation the user would enter the parameters as they would for the discrete step cellular automata including the probabilities for the birth and death events. Probabilities are used so that all experiments can be run with the same parameters and the simulator will convert them to rates as required for the algorithm. The conversion (line A2 in Figure \ref{fig:gill}) is done using the formula \[r = \lambda{e^{-\lambda{t}}}\]
where $\lambda$ is the probability of the event occurring and t is the time step (in this case 1). This is then multiplied by the population size for that species and all the rates are summed together (line A3 in Figure \ref{fig:gill}) and used as the total rate $r_{T}$ in the formula shown below. This is to maintain a memoryless process.

The time ($\tau$) until the next event is to occur (line C1 in Figure \ref{fig:gill}) is determined using the formula
\[\tau = \frac{-\ln{(1 - \text{Random Double})}}{r_{T}}\]

Determining the rates for each event is done in section A of Figure \ref{fig:gill} which is then used to determine the probability of each event occurring (section B). This is the equivalent of placing each event proportionally over the line in Figure \ref{fig:numline}. The event to occur is then chosen by generating a random number (line C3) and then determining which event it relates to (section D). Once the event is chosen, an organism for it to be applied to must be selected. This is done by randomly indexing the grid until a cell containing an organism of the correct species is found and the event is then executed on that organism. The time $\tau$ is also added to the total time (line C2 in Figure \ref{fig:gill}) and this process is repeated until the overall time crosses over into a new whole number at which point a snapshot of the system is taken and is counted as one step. For example, if T = 5.96 at the end of one event and then the value $\tau$ was chosen to be 0.05 (so T = 6.1), at the end of that event the state of the system would be represented to the user as one step. This would then continue until an event occurred which made T $>$ 7.0, then the new state of the system (current population sizes, etc.) would then be displayed to the user.

\subsection{Spatial Indices}
\label{sec:meth.impl.spat}

\subsubsection*{Patches}
\label{sec:meth.impl.spat.patches}

A ``patch" is defined as ``a continuous region of occupied sites assuming as contiguous the eight closest sites", i.e. is a connected component of a graph where each cell in the grid is a vertex and an edge is formed between cells containing the same species within their own Moore neighbourhood. The number of organisms in the patch is defined as its size or index (see Figure \ref{fig:patchEx}).

\begin{figure}[H]
\begin{center}
\includegraphics[scale = 0.91]{images/patchEx35.png}
\end{center}
\caption{Organisms in cells B2, C3 and D2 form one patch of size 3 and B5 is a patch of 1}
\label{fig:patchEx}
\end{figure}

The algorithm for measuring the number of patches and the biggest patch size is shown in Figure \ref{fig:pat}. In order to measure the number of patches without evaluating any cell twice, a check is kept whether the cell has been previously evaluated (line B4 in Figure \ref{fig:pat}). Using this method, the size of the patch is found using a breadth first search. As this process is run, the size of the largest patch is also maintained. For example, the process used to determine the biggest patch size and the number of patches for the grid in Figure \ref{fig:patchEx} can be found in Figure \ref{fig:patchExample}.

\begin{figure}[H]
\centering \small
\begin{footnotesize}
\begin{tabular}{ll}
A1 & \textbf{for} cells in grid\\
A2 & \hspace{0.25cm} \textbf{if} cell is not empty \textbf{and} cell has not been previously checked\\
A3 & \hspace{0.5cm} add cell to list if it is not already there\\
A4 & \hspace{0.5cm} noOfPatches = noOfPatches + 1\\
B1 & \hspace{0.75cm} \textbf{while} list not empty\\
B2 & \hspace{1.0cm} remove first cell from list\\
B4 & \hspace{1.0cm} mark cell as checked\\
B4 & \hspace{1.0cm} add previously unchecked neighbours containing same species to list\\
B5 & \hspace{1.0cm} patchSize = patchSize + 1\\
C1 & \hspace{0.75cm} \textbf{if} patchSize \textbf{$>$} biggestPatch\\
C2 & \hspace{1.0cm} biggestPatch = patchSize\\
\end{tabular}
\end{footnotesize}
\caption{General pseudocode for the algorithm for measuring patches}
\label{fig:pat}
\end{figure}

\begin{figure}[H]
\centering \small
\begin{footnotesize}
\begin{tabular}{lll}
1 & for cells A1 - B1, nothing happens. & \textbf{list = [], checked = []}\\ \\
2 & for cell B2 an unchecked, non-empty cell is & \textbf{list = [B2], checked = []}\\
 &  found so added to list. & \\ \\
2.1 & remove first item in list and mark as checked. & \textbf{list = [], checked = [B2]}\\
 & & \textbf{biggestPatch = 0, noOfPatches = 1, currentPatch = 1}\\ \\
2.2 & add cells B3 and C3 to list. & \textbf{list = [B3, C3], checked = [B2]}\\
 & & \textbf{biggestPatch = 0, noOfPatches = 1, currentPatch = 1}\\ \\
2.3 & remove first item in list and mark as checked. & \textbf{list = [C3], checked = [B2, B3]}\\
 & & \textbf{biggestPatch = 0, noOfPatches = 1, currentPatch = 2}\\ \\
2.4 & no new cells to add to list. & \textbf{list = [C3], checked = [B2, B3]}\\
 & & \textbf{biggestPatch = 0, noOfPatches = 1, currentPatch = 2}\\ \\
2.5 & remove first item in list and mark as checked. & \textbf{list = [], checked = [B2, B3, C3]}\\
 & & \textbf{biggestPatch = 0, noOfPatches = 1, currentPatch = 3}\\ \\
2.6 & no new cells to add to list. & \textbf{list = [], checked = [B2, B3, C3]}\\
 & & \textbf{biggestPatch = 0, noOfPatches = 1, currentPatch = 3}\\ \\
2.7 & biggestPatch = patchSize & \textbf{list = [], checked = [B2, B3, C3]}\\
 & & \textbf{biggestPatch = 3, noOfPatches = 1, currentPatch = 3}\\ \\
3 & cell B3 contains a checked cell so is skipped & \textbf{list = [], checked = [B2, B3, C3]}\\
 &  and cell B4 is empty so nothing happens. & \textbf{biggestPatch = 3, noOfPatches = 1, currentPatch = 3}\\ \\
4 & for cell B5 an unchecked, non-empty cell is & \textbf{list = [B5], checked = [B2, B3, C3]}\\
 & found so added to the list. & \textbf{biggestPatch = 3, noOfPatches = 1, currentPatch = 3}\\ \\
4.1 & remove first item in list and mark as checked. & \textbf{list = [], checked = [B2, B3, C3, B5]}\\
 & & \textbf{biggestPatch = 3, noOfPatches = 2, currentPatch = 1}\\ \\
4.2 & no new cells to add to list. & \textbf{list = [], checked = [B2, B3, C3, B5]}\\
 & & \textbf{biggestPatch = 3, noOfPatches = 2, currentPatch = 1}\\ \\
5 & for cells C1 - E5 nothing happens. & \textbf{list = [], checked = [B2, B3, C3, B5]}\\
 & & \textbf{biggestPatch = 3, noOfPatches = 2, currentPatch = 1}\\
6 & Returns biggestPatch as 3 and noOfPatches as 2 & \\
\end{tabular}
\end{footnotesize}
\caption{Worked example for determining the largest patch size and number of patches for the grid in Figure \ref{fig:patchEx}}
\label{fig:patchExample}
\end{figure}

\subsubsection*{Moran's I}
\label{sec:meth.impl.spat.moran}

Moran's I is a spatial autocorrelation index which is used to determine the spatial behaviour of the species of the system in regards to how ``clustered" they are. It looks to see how likely it is a that a cell is surrounded by other cells of the same type (i.e. if the current cell is occupied by an organism of species 1, how likely is it that the neighbouring cells are also occupied with organisms of species 1). This is done using the formula:
\[
I = \frac{N}{\sum_{i}\sum_{j} w_{ij}} \frac{\sum_{i} \sum_{j}(X_{i} - \bar{X}) (X_{j} - \bar{X})}{\sum_{i}(X_{i} - \bar{X})^2}
\]\

where:
\begin{itemize}
\item \emph{N} is the number of cells in the grid.
\item \emph{i} and \emph{j} are indexes of the grid.
\item \emph{$w_{ij}$} is the connection between cell \emph{i} and \emph{j}.
\item \emph{$X_{i}$} takes the value of 1 if the cell \emph{i} is occupied, or 0 otherwise.
\item \emph{$X_{j}$} takes the value of 1 if the cell \emph{i} and \emph{j} are of the same type, or 0 otherwise.
\item \emph{$\bar{X}$} = $\frac{\text{population size}}{\text{size of the grid}}$
\end{itemize}

The equation
\[
\frac{N}{\sum_{i}\sum_{j} w_{ij}}
\]

is the population divided by the total number of connections there are for each cell (E2 in Figure \ref{fig:moran}) i.e.:
\[
\frac{N}{8N}
\]

if the Moore neighbourhood is used. The formula
\[
\sum_{i} \sum_{j}(X_{i} - \bar{X}) (X_{j} - \bar{X})
\]
determines the spatial autocorrelation of each of the connections. This is done by summing up the values for every neighbouring of whether the cells are the same type or not (section C in Figure \ref{fig:moran}). This is then all divided by the total variation which is calculated by summing the value of each cell (1 if occupied, 0 otherwise) minus the mean value multiplied by itself (section D in Figure \ref{fig:moran}). This calculation returns a number between -1 (which indicates that every cell is neighboured by a cell which is of a different type) and a value close to 1 (which indicates that every cell is neighboured by one which is of the same type as it). A fully mixed system would give a value close to 0.

\begin{figure}[H]
\centering \small
\begin{footnotesize}
\begin{tabular}{ll}
A1 & \textbf{for} cell in grid:\\
A2 & \hspace{0.25in} xBar = xBar + cell value\\
A3 & xBar = xBar / grid length\\
B1 & \textbf{for} cell in grid:\\
B2 & \hspace{0.25in} \textbf{for} neighbour in neighbourhood:\\
C1 & \hspace{0.5in} \textbf{if} cell \textbf{and} neighbour are empty:\\
C2 & \hspace{0.75in} sac = sac + ((0 - xBar) * (0 - xBar))\\
C3 & \hspace{0.5in} \textbf{else if} cell \textbf{==} neighbour:\\
C4 & \hspace{0.75in} sac = sac + ((1 - xBar) * (1 - xBar))\\
C5 & \hspace{0.5in} \textbf{else}:\\
C6 & \hspace{0.75in} sac = sac + ((1 - xBar) * (0 - xBar))\\
D1 & \hspace{0.25in} \textbf{if} cell is empty:\\
D2 & \hspace{0.5in} tv = tv + ((0 - xBar) * (0 - xBar))\\
D3 & \hspace{0.25in} \textbf{else}:\\
D4 & \hspace{0.5in} tv = tv + ((1 - xBar) * (0 - xBar))\\
E1 & n = grid length\\
E2 & connections = n * size of neighbourhood\\
E3 & moransI = (n/connections) * (sac/tv)\\
\end{tabular}
\end{footnotesize}
\caption{Pseudocode for calculating Moran's I}
\label{fig:moran}
\end{figure}

\section{GUI}
\label{sec:meth.GUI}
The Graphical User Interface is split into two parts, the Grid and the Chart. Each part is optionally shown using the flags in the command line arguments. The flags are `-g' and `-c' for the grid and the chart respectively.

\subsection{Grid}
\label{sec:meth.GUI.grid}
The grid plots the position of each organism in the array so that the user can examine the spatial organisation of each of the species over time. For the grid, I used a window of a set size and set the dimensions of each cell based upon the size of the grid chosen by the user to maximise usage of the window. Each species within the simulation is given a different colour and each cell with an organism occupying it is shaded in the colour of that species or is otherwise left white (see Figure \ref{fig:GUI_Grid}). The grid is updated after each step by running through a stored list of changes (storing which cell has changed and what it has changed to) from the previous step and then changes the colours of the affected cells. For the implementation of the grid I used the drawing package `StdDraw' written by Robert Sedgewick (available at \url{http://introcs.cs.princeton.edu/java/stdlib/StdDraw.java.html}).

\begin{figure}[h]
\begin{center}
\includegraphics[scale = 0.75]{images/GUI_Grid.png}
\end{center}
\caption{The system with two species}
\label{fig:GUI_Grid}
\end{figure}

\subsection{Chart}
\label{sec:meth.GUI.chart}
The chart was developed to show the state of the system in real time (as it runs) and shows population size(see Figure \ref{fig:GUI_Chart}). It plots the population size of each species on a line chart (the colour of the line for each species matching the colour used in the grid) which updates at every step. The axes of the chart are re-scaled at each step so that the chart utilises the entire window, allowing as much detail to be shown as possible. This way the population sizes of each species (the y-axis) are shown relative to each other (the maximum shown is equivalent to the highest peak in population of any of the species) as opposed to relative to the size of the whole grid. This is useful for when any, or all, of the species only take up a small proportion of the entire grid, otherwise a lot of detail could be lost. For the x-axis(time), this is useful as the chart shows the state of the system as it is in that moment with how it got to that point, rather than losing detail to show part of the grid as empty for the steps which are yet to occur or may not occur at all (in the event of system extinction). 

\begin{figure}[H]
\begin{center}
\includegraphics[scale = 0.75]{images/GUI_Chart.png}
\end{center}
\caption{A snapshot of the chart with two species}
\label{fig:GUI_Chart}
\end{figure}

\subsection*{Summary}
This chapter showed all of the features of the simulator as well as how to use it. It also explained how the cellular automata and the Gillespie simulator were coded such that the reader should be able to recreate their own version. An explanation of the different metrics implemented for the experiments were also introduced and explained here, as was the pseudocode and some worked examples so that the reader would be able to understand how they are calculated.
